SUBROUTINE DG_MOT_PHI
!-----------------------------------------------------------------------
!
! Definition von G:  Drehimpulsgleichung
!

      use primvar,  only : DG, X, XZ, XA, MR, MD, MUr, MUp, ME, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use physco,   only : z1, z2, z12, pi, sqrt2pi, grav, clight
      use config,   only : np, BCflag
      use global,   only : tst, zz
      use geomvar,  only : S_flux, S_volZ, S_vol, S_volA
      use advecvar, only : momp_adv, dmomp_adv_dim2, dmomp_adv_dim1, dmomp_adv_di, dmomp_adv_dip1
      use matvar,   only : Pgas0Z, dPgas0ZdE, dPgas0ZdD, OPArosZ, dOPArosZdE, dOPArosZdD
      use viscvar,  only : divU, divU_flag, alpha_p, q_quad_p, muQ_p, Uterm_p
      use RBvar,    only : M_star
      use zvar,     only : H_pZ, dH_pZdE, dH_pZdR, H_p, dH_pdE, dH_pdR


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      DG(MUp,MUp,DG_i,1)      =  z1
      DG(MUp,MR,DG_i,1)       =  z12 * sqrt(grav*M_star/X(MR,1)**3)
      DG(MUp,MUp,DG_i,2)      =  z1
      DG(MUp,MR,DG_i,2)       =  z12 * sqrt(grav*M_star/X(MR,2)**3)

   ! aeussere Pseudozellen: i=np, i=np-1 & i=np-2
      DG(MUp,MUp,DG_i,np)     =  z1
!~       DG(MUp,MR,DG_i,np)      =  z12 * sqrt(grav*M_star/X(MR,np)**3)
      DG(MUp,MUp,DG_i,np-1)   =  z1
!~       DG(MUp,MR,DG_i,np-1)    =  z12 * sqrt(grav*M_star/X(MR,np-1)**3)
      DG(MUp,MUp,DG_i,np-2)   =  z1
!~       DG(MUp,MR,DG_i,np-2)    =  z12 * sqrt(grav*M_star/X(MR,np-2)**3)


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

DG(MUp,MR,DG_im1,i)=(muQ_p(i)*pi*XZ(MR,i)**2*tst*(z1  - BCflag(i) - BCflag(i-1))*(XZ(MUp,i) - XZ(MUp,i-1))*z12 &
*z2 *zz ) / (XZ(MR,i+1) - XZ(MR,i-1))**2
DG(MUp,MR,DG_i,i)=-(pi*X(MD,i)*(X(MR,i) + X(MR,i+1))**2*X(MUp,i)*z12 ) + X(MD,i)*S_vol(i)*X(MUp,i)*z12 *z2  - &
S_flux(i)*zz *momp_adv(i) - XZ(MR,i)*(-(pi*(XA(MR,i) + X(MR,i))) - pi*(X(MR,i) - XA(MR,i)) + pi*tst*XZ(MUr,i)*z2 *zz &
)*momp_adv(i) - pi*tst*z12 *(-(muQ_p(i)*(XZ(MUp,i) + XZ(MUp,i-1))*(z1  - BCflag(i) - BCflag(i-1))*z12 *zz ) - &
muQ_p(i)*XZ(MR,i)*Uterm_p(i)*z2 *zz  + (muQ_p(i+1)*XZ(MR,i+1)**2*(z1  - BCflag(i) - BCflag(i+1))*(XZ(MUp,i+1) - XZ(MUp,i))*z2 *zz &
) / (XZ(MR,i+2) - XZ(MR,i))**2 - XZ(MR,i)**2*Uterm_p(i)*(-(divU_flag(i)*divU(i)*q_quad_p**2*XZ(MD,i)*XZ(MR,i)*z2 *zz ) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i)*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(-(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz ) + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2 + &
(alpha_p(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_p(i)*XZ(MD,i)**2*z12 &
*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) - (alpha_p(i)*X(MD,i)*XZ(MD,i)*z12 &
*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)**2*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))))
DG(MUp,MR,DG_ip1,i)=pi*X(MD,i)*(X(MR,i) + X(MR,i+1))**2*X(MUp,i)*z12  + X(MD,i)*S_vol(i)*X(MUp,i)*z12 *z2  + &
S_flux(i+1)*zz *momp_adv(i+1) + XZ(MR,i+1)*(-(pi*(XA(MR,i+1) + X(MR,i+1))) - pi*(X(MR,i+1) - XA(MR,i+1)) + pi*tst*XZ(MUr,i+1)*z2 &
*zz )*momp_adv(i+1) - pi*tst*z12 *(muQ_p(i+1)*(XZ(MUp,i) + XZ(MUp,i+1))*(z1  - BCflag(i) - BCflag(i+1))*z12 *zz  + &
(muQ_p(i)*XZ(MR,i)**2*(z1  - BCflag(i) - BCflag(i-1))*(XZ(MUp,i) - XZ(MUp,i-1))*z2 *zz ) / (XZ(MR,i+1) - XZ(MR,i-1))**2 + &
muQ_p(i+1)*XZ(MR,i+1)*(z1  - BCflag(i) - BCflag(i+1))*z2 *(-(((XZ(MUp,i) + XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - &
XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - XZ(MR,i)))*zz  - &
XZ(MR,i)**2*Uterm_p(i)*(-((divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i+1)*z2 *zz ) / (BCflag(i) + S_volZ(i))) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz  + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2) + XZ(MR,i+1)**2*(z1  - BCflag(i) - &
BCflag(i+1))*(-(((XZ(MUp,i) + XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - &
XZ(MR,i)))*(-(divU_flag(i+1)*divU(i+1)*q_quad_p**2*XZ(MD,i+1)*XZ(MR,i+1)*z2 *zz ) + &
(divU_flag(i+1)*pi*q_quad_p**2*XZ(MD,i+1)*XZ(MR,i+1)**2*XZ(MUr,i+1)*z2 *zz ) / (BCflag(i+1) + S_volZ(i+1)) + &
(alpha_p(i+1)*XZ(MD,i+1)*Sqrt(Pgas0Z(i+1))*dH_pZdR(i+1)) / Sqrt(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1))) + &
(alpha_p(i+1)*XZ(MD,i+1)**2*z12 *Sqrt(Pgas0Z(i+1))*dH_pZdR(i+1)) / (sqrt2pi*(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1)))**1.5*H_pZ(i+1)) - &
(alpha_p(i+1)*X(MD,i+1)*XZ(MD,i+1)*z12 *H_pZ(i+1)*dH_pdR(i+1)*dPgas0ZdD(i+1)) / (sqrt2pi*H_p(i+1)**2*Sqrt(XZ(MD,i+1) / &
(sqrt2pi*H_pZ(i+1)))*Sqrt(Pgas0Z(i+1)))))
DG(MUp,MR,DG_ip2,i)=-(pi*tst*z12 *(-((muQ_p(i+1)*XZ(MR,i+1)**2*(z1  - BCflag(i) - BCflag(i+1))*(XZ(MUp,i+1) - &
XZ(MUp,i))*z2 *zz ) / (XZ(MR,i+2) - XZ(MR,i))**2) - (divU_flag(i+1)*pi*q_quad_p**2*XZ(MD,i+1)*XZ(MR,i+1)**4*XZ(MUr,i+2)*(z1  - &
BCflag(i) - BCflag(i+1))*z2 *(-(((XZ(MUp,i) + XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - &
XZ(MR,i)))*zz ) / (BCflag(i+1) + S_volZ(i+1))))
DG(MUp,MD,DG_im2,i)=-(XZ(MR,i)*S_flux(i)*XZ(MUp,i-2)*zz *dmomp_adv_dim2(i))
DG(MUp,MD,DG_im1,i)=-(XZ(MR,i)*S_flux(i)*XZ(MUp,i-1)*zz *dmomp_adv_dim1(i)) + &
XZ(MR,i+1)*S_flux(i+1)*XZ(MUp,i-1)*zz *dmomp_adv_dim2(i+1)
DG(MUp,MD,DG_i,i)=pi*(X(MR,i) + X(MR,i+1))**2*X(MUp,i)*(X(MR,i+1) - X(MR,i))*z12  + &
pi*XZ(MR,i)**2*tst*Uterm_p(i)*z12 *(-(divU_flag(i)*divU(i)*q_quad_p**2*XZ(MR,i)**2*zz ) - (alpha_p(i)*XZ(MD,i)*z12 *zz &
*Sqrt(Pgas0Z(i))) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5) + (alpha_p(i)*zz *H_pZ(i)*Sqrt(Pgas0Z(i))) / Sqrt(XZ(MD,i) / &
(sqrt2pi*H_pZ(i))) + (alpha_p(i)*XZ(MD,i)*z12 *H_pZ(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)*Sqrt(XZ(MD,i) / &
(sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))) - XZ(MR,i)*S_flux(i)*XZ(MUp,i)*zz *dmomp_adv_di(i) + XZ(MR,i+1)*S_flux(i+1)*XZ(MUp,i)*zz &
*dmomp_adv_dim1(i+1)
DG(MUp,MD,DG_ip1,i)=-(pi*XZ(MR,i+1)**2*tst*(z1  - BCflag(i) - BCflag(i+1))*z12 *(-(((XZ(MUp,i) + &
XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - &
XZ(MR,i)))*(-(divU_flag(i+1)*divU(i+1)*q_quad_p**2*XZ(MR,i+1)**2*zz ) - (alpha_p(i+1)*XZ(MD,i+1)*z12 *zz *Sqrt(Pgas0Z(i+1))) / &
(sqrt2pi*(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1)))**1.5) + (alpha_p(i+1)*zz *H_pZ(i+1)*Sqrt(Pgas0Z(i+1))) / Sqrt(XZ(MD,i+1) / &
(sqrt2pi*H_pZ(i+1))) + (alpha_p(i+1)*XZ(MD,i+1)*z12 *H_pZ(i+1)*dPgas0ZdD(i+1)) / (sqrt2pi*H_p(i+1)*Sqrt(XZ(MD,i+1) / &
(sqrt2pi*H_pZ(i+1)))*Sqrt(Pgas0Z(i+1))))) - XZ(MR,i)*S_flux(i)*XZ(MUp,i+1)*zz *dmomp_adv_dip1(i) + &
XZ(MR,i+1)*S_flux(i+1)*XZ(MUp,i+1)*zz *dmomp_adv_di(i+1)
DG(MUp,MD,DG_ip2,i)=XZ(MR,i+1)*S_flux(i+1)*XZ(MUp,i+2)*zz *dmomp_adv_dip1(i+1)
DG(MUp,ME,DG_i,i)=pi*XZ(MR,i)**2*tst*Uterm_p(i)*z12 *((alpha_p(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdE(i)) / &
Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_p(i)*XZ(MD,i)**2*z12 *Sqrt(Pgas0Z(i))*dH_pZdE(i)) / (sqrt2pi*(XZ(MD,i) / &
(sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) + (alpha_p(i)*XZ(MD,i)*z12 *H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + &
dPgas0ZdE(i))) / (Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))
DG(MUp,ME,DG_ip1,i)=-(pi*XZ(MR,i+1)**2*tst*(z1  - BCflag(i) - BCflag(i+1))*z12 *(-(((XZ(MUp,i) + &
XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - &
XZ(MR,i)))*((alpha_p(i+1)*XZ(MD,i+1)*Sqrt(Pgas0Z(i+1))*dH_pZdE(i+1)) / Sqrt(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1))) + &
(alpha_p(i+1)*XZ(MD,i+1)**2*z12 *Sqrt(Pgas0Z(i+1))*dH_pZdE(i+1)) / (sqrt2pi*(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1)))**1.5*H_pZ(i+1)) + &
(alpha_p(i+1)*XZ(MD,i+1)*z12 *H_pZ(i+1)*(-((X(MD,i+1)*dPgas0ZdD(i+1)*dH_pdE(i+1)) / (sqrt2pi*H_p(i+1)**2)) + dPgas0ZdE(i+1))) / &
(Sqrt(XZ(MD,i+1) / (sqrt2pi*H_pZ(i+1)))*Sqrt(Pgas0Z(i+1)))))
DG(MUp,MUr,DG_i,i)=(divU_flag(i)*pi**2*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**5*tst*Uterm_p(i)*z12 *z2 *zz ) / &
(BCflag(i) + S_volZ(i)) - pi*XZ(MR,i)**2*tst*z2 *zz *momp_adv(i)
DG(MUp,MUr,DG_ip1,i)=-(pi*tst*z12 *((divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**4*XZ(MR,i+1)*Uterm_p(i)*z2 &
*zz ) / (BCflag(i) + S_volZ(i)) + (divU_flag(i+1)*pi*q_quad_p**2*XZ(MD,i+1)*XZ(MR,i+1)**5*(z1  - BCflag(i) - BCflag(i+1))*z2 &
*(-(((XZ(MUp,i) + XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - XZ(MR,i)))*zz ) / &
(BCflag(i+1) + S_volZ(i+1)))) + pi*XZ(MR,i+1)**2*tst*z2 *zz *momp_adv(i+1)
DG(MUp,MUr,DG_ip2,i)=(divU_flag(i+1)*pi**2*q_quad_p**2*XZ(MD,i+1)*XZ(MR,i+1)**4*XZ(MR,i+2)*tst*(z1  - BCflag(i) &
- BCflag(i+1))*z12 *z2 *(-(((XZ(MUp,i) + XZ(MUp,i+1))*z12 ) / XZ(MR,i+1)) + ((XZ(MUp,i+1) - XZ(MUp,i))*z2 ) / (XZ(MR,i+2) - &
XZ(MR,i)))*zz ) / (BCflag(i+1) + S_volZ(i+1))
DG(MUp,MUp,DG_im2,i)=-(XZ(MD,i-2)*XZ(MR,i)*S_flux(i)*zz *dmomp_adv_dim2(i))
DG(MUp,MUp,DG_im1,i)=muQ_p(i)*pi*XZ(MR,i)**2*tst*(z1  - BCflag(i) - BCflag(i-1))*z12 *(-((z12 *zz ) / XZ(MR,i)) &
- (z2 *zz ) / (XZ(MR,i+1) - XZ(MR,i-1))) - XZ(MD,i-1)*XZ(MR,i)*S_flux(i)*zz *dmomp_adv_dim1(i) + &
XZ(MD,i-1)*XZ(MR,i+1)*S_flux(i+1)*zz *dmomp_adv_dim2(i+1)
DG(MUp,MUp,DG_i,i)=pi*X(MD,i)*(X(MR,i) + X(MR,i+1))**2*(X(MR,i+1) - X(MR,i))*z12  - pi*tst*z12 &
*(muQ_p(i+1)*XZ(MR,i+1)**2*(z1  - BCflag(i) - BCflag(i+1))*(-((z12 *zz ) / XZ(MR,i+1)) - (z2 *zz ) / (XZ(MR,i+2) - XZ(MR,i))) - &
muQ_p(i)*XZ(MR,i)**2*(z1  - BCflag(i) - BCflag(i-1))*(-((z12 *zz ) / XZ(MR,i)) + (z2 *zz ) / (XZ(MR,i+1) - XZ(MR,i-1)))) - &
XZ(MD,i)*XZ(MR,i)*S_flux(i)*zz *dmomp_adv_di(i) + XZ(MD,i)*XZ(MR,i+1)*S_flux(i+1)*zz *dmomp_adv_dim1(i+1)
DG(MUp,MUp,DG_ip1,i)=-(muQ_p(i+1)*pi*XZ(MR,i+1)**2*tst*(z1  - BCflag(i) - BCflag(i+1))*z12 *(-((z12 *zz ) / &
XZ(MR,i+1)) + (z2 *zz ) / (XZ(MR,i+2) - XZ(MR,i)))) - XZ(MD,i+1)*XZ(MR,i)*S_flux(i)*zz *dmomp_adv_dip1(i) + &
XZ(MD,i+1)*XZ(MR,i+1)*S_flux(i+1)*zz *dmomp_adv_di(i+1)
DG(MUp,MUp,DG_ip2,i)=XZ(MD,i+2)*XZ(MR,i+1)*S_flux(i+1)*zz *dmomp_adv_dip1(i+1)

      end do


END SUBROUTINE DG_MOT_PHI
